Sparsity-based method for ring artifact elimination in computed tomography

Ring artifact elimination is one of the popular problems in computed tomography (CT). It appears in the reconstructed image in the form of bright or dark patterns of concentric circles. In this paper, based on the compressed sensing theory, we propose a method for eliminating the ring artifact during the image reconstruction. The proposed method is based on representing the projection data by a sum of two components. The first component contains ideal correct values, while the latter contains imperfect error values causing the ring artifact. We propose to minimize some sparsity-induced norms corresponding to the imperfect error components to effectively eliminate the ring artifact. In particular, we investigate the effect of using different sparse models, i.e. different sparsity-induced norms, on the accuracy of the ring artifact correction. The proposed cost function is optimized using an iterative algorithm derived from the alternative direction method of multipliers. Moreover, we propose improved versions of the proposed algorithms by incorporating a smoothing penalty function into the cost function. We also introduce angular constrained forms of the proposed algorithms by considering a special case as follows. The imperfect error values are constant over all the projection angles, as in the case where the source of ring artifact is the non-uniform sensitivity of the detector. Real data and simulation studies were performed to evaluate the proposed algorithms. Results demonstrate that the proposed algorithms with incorporating smoothing penalty and their angular constrained forms are effective in ring artifact elimination.


Introduction
Ring artifact is one of the well-known image degradations that restrict the quality of X-ray computed tomography (CT) images [1]. In addition to general-purpose medical CT scanners, it also occurs in CT images that are generated from a flat-panel detector such as C-arm CT, dental CT, and micro-CT. The ring artifact can be recognized in the image domain as bright or dark patterns of concentric circles. These circles vary in intensity value from ring to ring. Also, it can be interpreted in the sinogram domain as bright or dark vertical lines having a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 different intensity values. The ring artifact mainly arises as a consequence of some physical factors associated with non-uniform detector sensitivity, such as 1) existence of imperfect or dead detector elements, 2) incomplete calibration of the detector elements, 3) impurity or dust in the scintillator screen, or 4) presence of variation in the hardness or intensity of x-ray beams. The ring artifact may resemble some structure of the scanned object, which prevents the resulting image from being used in various applications. Hence, the elimination of ring artifact is a significant issue. One conventional technique that is used in ring artifact elimination is called flat-field correction. It was reported that such simple technique could eliminate ring artifact partially [2]. Other methods that deal with the ring artifact problem can be divided into hardware-and software-based methods. The removal of ring artifact in most software-based methods was performed before or after the reconstruction process. However, such methods may introduce other artifacts. On the other hand, the elimination of the ring artifact during the reconstruction overcomes the limitations of such methods. In this paper, inspired by the compressed sensing (CS) theory [3], we propose a software-based method to estimate and eliminate the ring artifact during the image reconstruction process. Some sparsity-based methods were presented to suppress the ring artifacts [4][5][6][7][8]. In this work, we introduce using sparsityinduced norms like ℓ 0 , Huber ℓ 0 , and Huber ℓ 1 in ring artifact removal through the reconstruction steps. The sinogram data corresponding to the image containing the ring artifact is also contaminated by incorrect measurement values similar to the example shown in Fig 1. The sinogram data in such case is composed of the perfect ideal values and the imperfect error values causing the ring artifact. The contribution of this paper is as follows. First, we propose to use some sparsity-induced norms of the imperfect component as data fidelity to enforce the sparsity of the imperfect component to eliminate the ring artifact. Second, we investigate the effect of using different sparsity-induced norms on the ring artifact elimination. Third, we propose to improve image quality in the proposed algorithms by incorporating a smoothing penalty function into the cost function. Fourth, we study a particular case of the developed algorithms, where the values of the imperfect error component are constant over all the projection angles. Finally, we remark that this work can be considered a particular application of the use of sparsity-induced norms in data fidelity to eliminate the effect of abnormal error data in the sinogram [9,10]. Such approach was proposed in our previous work together with a special application to the metal artifact reduction problem [11]. This paper demonstrates that the similar approach works well for the ring artifact correction problem. This paper is organized as follows. The related work is presented in Section 2. The proposed algorithms are detailed in Section 3. After that, we show the results of the experimental studies in Section 4. Finally, the conclusion is provided in 5.

Related work
There are numerous methods that have been investigated to suppress the ring artifact. These methods can be classified into two general categories. The first category contains methods based on the modification of hardware, such as those detailed in [12,13]. The second category contains software-based methods. Due to the necessity of modifying the device hardware and consequently the extra cost in the first class, most of the investigated methods belong to the second category. The software-based methods can be further divided into three classes as follows: • Pre-processing methods: remove the vertical stripes from the sinogram before image reconstruction.
• Post-processing methods: remove the ring artifact from the image after image reconstruction.
• In-processing methods: remove the ring artifact during image reconstruction.

Pre-processing methods
A technique based on a combination of wavelet decomposition and Fourier filtering was investigated to eliminate the stripes from the sinogram [14]. Also, a sinogram-based method was proposed by supposing smoothness of the detector measurements for every projection angle [15]. A derivative-based algorithm was developed to detect the mis-calibrated and defective detector elements [16]. This algorithm used a bias correction model and a weighted moving average filter to remove the stripe artifacts. An approach focused on filtering scheme by morphological operations was developed for the ring artifact elimination [17]. A detector line ratio method based on computing the ratio between the neighboring detector elements was proposed in the sinogram space [18]. Furthermore, based on classifying the stripe artifacts in the sinogram into different types, three methods were designed to deal with each type individually [19]. Also, a combination of these three methods was proposed, and they concluded that this combination could suppress all types of artifacts. A ring artifact removal algorithm was proposed in the form of a 1D filter and can be combined with the standard filter used in the FBP method [20]. Using wavelet filtering and weighted moving average filter in the sinogram space, an approach was presented to remove the ring artifact [21]. The removal of stripe artifacts in the sinogram was performed via variable window averaging and weighted moving average filters [22]. A U-Net-like deep neural network was employed to reduce the stripe artifacts in the sinogram space [23]. However, filtering the sinogram to remove the vertical stripes may introduce new artifacts.

Post-processing methods
A correction algorithm was proposed by utilizing mean and median filters in both Cartesian and polar coordinate systems. They concluded that implementing their algorithm in the polar coordinate system gives higher performance [24]. A method based on generative adversarial networks was investigated to remove the ring artifact from the image in the Cartesian coordinates [25]. Moreover, several methods have been investigated to suppress the ring artifact from the image in the polar coordinate system. For example, an iterative method was developed to remove the ring artifact from CT images utilizing relative total variation (TV) [26]. A method using independent component analysis was presented for suppressing the ring artifact [27]. This method divides the image into many independent components. Then, the components that include the rings are identified to apply a smoothing filter only on the ring components. A polar wavelet Gaussian filtering technique was presented to suppress the ring artifact [28]. Furthermore, a method that depends on TV-Stokes and unidirectional TV model was introduced to eliminate the ring artifact from cone-beam CT images [29]. Based on a unidirectional relative variation model, an approach was investigated to suppress the ring artifact [30]. Also, a variation-based method was formulated using ℓ 0 -norm/ℓ 1 -norm regularization term to remove the ring artifact from the reconstructed images [4]. The reduction of ring artifact was achieved based on sparse domain regularized stripe decomposition technique and guided image filtering [5]. Moreover, a radial basis function neural network was designed to eliminate the ring artifact from the reconstructed images [6]. Using generative adversarial networks and unidirectional relative TV, a method was presented to suppress the artifact from the CT images [31]. Most of the post-processing methods consist of the following steps. Transform the image contaminated with the ring artifact to the polar coordinates and carry out the filtering, then return to the original Cartesian coordinates. The rationale behind using such transformation is as follows. The ring artifact in the Cartesian coordinates can be converted into lines in the polar coordinates, and filtering out the straight lines is much easier than the rings. However, image quality achieved by these methods relies on the accurate determination of the center of the ring artifact. This is because a small difference in the center determination could convert the rings into curvy lines instead of the straight lines. Such curvy lines distort the output image in the Cartesian coordinates.

In-processing methods
An optimization method based on TV and dictionary learning was investigated to correct the ring artifact [7]. Recently, an algorithm for the ring artifact reduction was proposed using ring TV and correction coefficient [8]. The proposed method follows this class of methods, where removing the ring artifact during the reconstruction process does not require any data modification before or after the reconstruction. Also, the probability of introducing other distortions like the pre-or post-processing methods is low.

Problem formulation
In the case of ring artifact existence in the image, the image reconstruction model in x-ray CT can be described as follows: Ax ¼b þl; ð1Þ Consequently, the imperfect vectorl can be considered as a sparse vector. The sparsity of imperfect vectorl inspires us to use the CS theory in ring artifact elimination problem. We exploit the CS theory in our work in order to estimate the imperfect data that cause the ring artifact and eliminate them from the projection data during image reconstruction to obtain a ring-free image. Therefore, we propose to formulate a cost function in image reconstruction by utilizing some sparsity-induced norms ofl as data fidelity to enforce the sparsity of imperfect component as follows: Minimize ðx;lÞ SðlÞ subject to Ax Àb ¼l; ð2Þ where SðlÞ can be one of well-known sparse models in CS such as ℓ 1 -norm or ℓ 0 -norm. The vectorl in Eq (2) is described to include only the non-ideal elements that generate the ring artifact. However, in real CT imaging, both the statistical noise and imperfect error values are included together inl. In such case, Huber loss function can be utilized instead of ℓ 1 -norm as a sparse model because it is a hybrid of ℓ 1 -norm and ℓ 2 -norm. Throughout the paper, we call it Huber ℓ 1 -norm. By the same manner, we suggest defining Huber ℓ 0 -norm as a combination of ℓ 0 -norm and ℓ 2 -norm. In summary, SðlÞ can be one of the following functions: SðlÞ ¼ klk SðlÞ ¼ SðlÞ ¼ Generally, it is clear from the different causes of ring artifact that the values of the imperfect vector in the projection data may be dependent on the measured angles. In other words, the imperfect values are variables for all the angles. In addition to this general case, we suggest studying a special case where the ring artifact originates from physical factors related to hardware, such as the existence of some imperfect detector bins. Then, the imperfect elements inl occur due to the imperfect detector bins. In CT, if some detector bins have errors, the incorrect values corresponding to those bins will be repeated for all the angles θ. Consequently, the measured value of the imperfect vectorl at any detector bin is constant for all the angles θ. Hereafter, such model having constant imperfect elements in the angular direction is called angular constrained form. It can be mathematically expressed as follows. Let us express the system matrix A using submatrices corresponding to each sinogram angle as where A i denotes the system matrix corresponding to the i-th projection angle. Similarly, we express the projection data asb ¼ ðb 1 ;b 2 ; . . . ;b Y Þ T , whereb i denotes the projection data corresponding to the i-th projection angle. Using these notations, the image reconstruction problem in the angular constrained form can be formulated as Minimize SðlÞ subject to A yx Àb y ¼l; 8y: ð7Þ We remark thatl is constant for all θ, which originates from the constraint that values of elements in the imperfect vector are same for all the angles.

Proposed algorithms
The formulated cost function can be optimized using any iterative algorithm. For example, first-order primal-dual method, alternative direction method of multipliers (ADMM) [32], and proximal splitting. In this work, we use ADMM method as a prototype example to develop iterative algorithms to solve the formulated cost function. We start with using ℓ 1 -norm as an example of the sparsity-induced norm SðlÞ in the proposed cost function as follows: Minimize ðx;lÞ klk 1 1 subject to Ax Àb ¼l ð8Þ The augmented Lagrangian corresponding to Eq (8) is whereũ denotes the Lagrange multiplier and ρ > 0 is a positive constant. Using the augmented Lagrangian function in Eq (9), ADMM method finds the solution to the primal-dual problem ðx;ũÞ according to the following three steps: The update of imagex is performed through the minimization in Eq (10). Since DðxÞ ¼ r 2 kl k À Ax kþ1 þb þũ k k 2 2 is a differentiable function in this case, it can be minimized using the standard gradient descent (GD) method as follows: where α > 0 denotes the step size, L is the maximum number of iterations in the GD method, and r x Dðx l Þ is calculated as follows: With respect to the update of the imperfect error vectorl, it can be done by the minimization in Eq (11). By recalling the definition of proximity operator as follows: the necessary computation in Eq (11) is reduced to Àb Àũ k and μ = 1/ρ. It can be solved directly using a simple formula in closed form via the soft thresholding [33] function as follows: where the soft-thresholding(�) is defined by The final iterative algorithm can be summarized in Algorithm 1. Algorithm 1 ℓ 1 -ring algorithm.
We followed the same steps to derive three iterative algorithms corresponding to the remaining three sparsity-induced norms mentioned above, i.e. ℓ 0 -norm, Huber ℓ 1 -norm and Huber ℓ 0 -norm. The result of the derivation shows the following. Almost same algorithms can be obtained for the cases of using four different norms, and the only difference among the four cases appears in the form of the proximity operator in Step 2. For clarity, the proximity operators corresponding to the four cases are summarized in Table 1. Throughout the paper, we call the algorithms corresponding to ℓ 1 -norm, ℓ 0 -norm, Huber ℓ 1 -norm, and Huber ℓ 0 -norm by ℓ 1 -ring, ℓ 0 -ring, Hℓ 1 -ring, and Hℓ 0 -ring algorithms, respectively.

Extension with smoothing penalty
We propose to incorporate a smoothing penalty HðxÞ into the proposed cost function in Eq (2) to strengthen the power to remove the artifact. The new cost function with the smoothing penalty is given by where β > 0 denotes the hyper-parameter that controls the strength of smoothing. In this work, we choose to use a smoothing penalty similar to Markov random fields (MRFs) [34]. This is because it possesses the advantage of edge-preserving property, and it is also a differentiable function that can be easily minimized. With respect to the optimization of penalized cost function, we first describe the ℓ 1 -norm case, followed by explaining how it can be changed to the other three norms. The cost function in the ℓ 1 -norm case is given by where HðxÞ is the smoothing penalty defined by where U j denotes the set of eight neighboring pixels around the j-th pixel, and w jj is the weighting parameter computed as the Euclidean distance between the j-th pixel and thej-th pixel. This modified cost function can be optimized by using iterative algorithms derived from ADMM similar to those in the previous section. The same ADMM steps yield, except Step 1, which can be modified and minimized as follows: x lþ1 ¼x l À arx Dðx l Þ À b 0 rx Hðx l Þ; l ¼ 1; 2; . . . ; L; where b 0 ¼ ab, and rx Hðx l Þ is calculated as follows: The final algorithm is detailed in Algorithm 2. Algorithm 2 ℓ 1 -smth-ring algorithm.
Input the projection datab. Set initial value forx 1 ,l 1 andũ 1 . Set the parameters (α > 0, ρ > 0 and b 0 ). for k = 1, 2, . . ., K dõ y 1 ¼x k for l = 1, 2, . . ., L dõ y lþ1 ¼ỹ l À arỹ Dðỹ l Þ À b 0 rỹ Hðỹ l Þ end x kþ1 ¼ỹ L l kþ1 ¼ Prox Similarly to the case without the smoothing penalty, the iterative algorithms corresponding to the other three norms mentioned above are given as follows: change the proximity operator in Algorithm 2, as shown in Table 1. Throughout the paper, we call the four obtained algorithms by ℓ 1 -smth-ring, ℓ 0 -smth-ring, Hℓ 1 -smth-ring, and Hℓ 0 -smth-ring algorithms, respectively. In addition, we call a collection of these four algorithms using the smoothing penalty by improved algorithms to distinguish them from the algorithms without the smoothing.

Constrained algorithms
This section describes iterative algorithms for the angular constrained form where elements of imperfect error vector have the same values for all the angles, i.e. the formulation of Eq (7). A simple algorithm for this special case can be given as follows. The proposed algorithms in the previous sections are only modified by including the constraint that the values of elements in the imperfect vectorl are constant for all the angles (the angular-independency constraint). We can observe the following with a deep look at Step 4 of the iterative algorithms at an iteration k. Each element λ rθ of the imperfect vectorl are updated for every detector bin (r) in every projection angle (θ) based on each elements d rθ of the vectord defined above as follows: for r = 1, 2, . . ., R do for θ = 1, 2, . . ., Θ do l kþ1 ry ¼ Prox Applying the angular-independency constraint means that the update ofl occurs only at every detector bin as follows: for r = 1, 2, . . ., R do l kþ1 r ¼ Prox Using the update ofl mentioned above in each iterative algorithm described in the previous sections, the corresponding iterative algorithm including the angular-independency can be obtained. We call the resulting algorithms by (cstr-ℓ 1 -ring, cstr-ℓ 0 -ring, cstr-Hℓ 1 -ring, and cstr-Hℓ 0 -ring algorithms) for the case without the smoothing, and (cstr-ℓ 1 -smth-ring, cstr-ℓ 0smth-ring, cstr-Hℓ 1 -smth-ring, and cstr-Hℓ 0 -smth-ring algorithms) for the case with the smoothing.

Convergence properties
The convergence of iterative algorithm using ℓ 1 -norm or Huber ℓ 1 -norm is guaranteed for the following reason. These algorithms can be considered as an application of ADMM to the convex minimization subject to the linear constraint. However, the iterative algorithms using ℓ 0norm or Huber ℓ 0 -norm are not guaranteed to converge due to the non-convex nature of cost functions. However, as we will show in the next section, the algorithms using ℓ 0 -norm or Huber ℓ 0 -norm converge to nice approximate solutions in practice. Also, their power to eliminate the ring artifact is much stronger compared to the ℓ 1 -norm based algorithms. This observation motivated us to keep the ℓ 0 -norm based algorithms in this paper in spite of the lack of mathematical guarantee in the convergence. We note that similar behaviors, i.e. approximate solution with the ℓ 0 -norm works better than the exact solution with the ℓ 1 -norm, have been observed in other application areas of CS.

Experimental setup
To evaluate the proposed algorithms, we performed experimental studies using real data in addition to simulation data. In the real data, some experiments for propagation-based phasecontrast imaging were performed by synchrotron radiation X-ray micro-CT (Spring-8, BL20XU beamline), Japan [35,36]. A sample of duplex stainless steel was scanned with an xray energy of 37.7 (KeV). The sample-to-detector distance was varied between 8 and 1200 (mm). Also, a CMOS camera of 4.0 megapixels, with a 10 (μm) thick scintillator was used for measuring the projection data. The number of projection data was 1800 over 180 degrees, and each projection data consisted of 1024 × 1024 pixels. From the conventional FBP reconstruction of this data, the central slice of the reconstructed 3D image was contaminated with strong ring artifact. Moreover, in the simulation data, we used a 2D slice obtained from a dataset of screening chest CT images. The dataset includes 68 volumes for various patients scanned using Hitachi CT-W950SR scanner. Each volume is composed of 25 to 31 transaxial slices, where each slice is composed of 320 × 320 (pixels) with a slice thickness of 10 (mm) and a pixel size of 1 × 1 (mm). The 2D chest slice was resized to 512 × 512 (pixels). The simulation studies were carried out with and without adding noise to the projection data. It is known that the statistical noise in the CT measurements follows Poisson distribution [37]. Assume that b 0 is the number of average photon counts transmitted from the x-ray source, and b i is the detected photons measured after attenuation by the scanned object at detector bin i. Then, the noisy projection data g i at detector bin i can be calculated using Eqs (24) and (25). Moreover, to simulate the ring artifact, we first converted the image into the sinogram by the forward projection. Then, the generated sinogram was altered by assuming that some selected detector bins have non-uniform sensitivity compared to normal bins. Finally, using the altered sinogram, we reconstructed the image by the standard filtered back-projection (FBP) method to obtain the corrupted image. For the purpose of performance evaluation, we compared our algorithms with the correction method [28] and the combination method [19]. We refer to these methods as M1 and M2, respectively. Three parameters are required to implement the M1 method (Gaussian filter width σ, decomposition level L, and wavelet base function W). On the other hand, to correct the image by the M2 method, some parameters need to be adjusted, such as the ratio between the defective and background values R. Also, the size of the smoothing filters that remove the small to medium and large stripes in the sinogram (S m , S l ) is required for executing the M2 method. In order to obtain satisfactory results with the M1 and M2 methods, we selected their parameters based on the best image quality in terms of RRMSE values besides visual evaluation. Furthermore, to study the effect of using the standard ℓ 2norm on ring artifact removal, we compared our proposed algorithms with two algorithms using ℓ 2 -norm for the data fidelity term with and without smoothing penalty. Additionally, the GD method was used for the optimization problem of these two algorithms, which we name ℓ 2 and ℓ 2 -smth algorithms. All the algorithms were implemented using C++ platform on a machine with Intel(R) core(TM) i7-2670QM CPU @ 2.20 GHz processor and 6 GB RAM. The average execution time of the proposed algorithms for every iteration was 6.0 seconds.

Image quality measurements
To evaluate image quality achieved by the proposed algorithms, we used two image quality metrics throughout the simulation studies. Suppose thatz ¼ ðz 1 ; . . . ; z N Þ > denotes the reference image andx denotes the reconstructed or corrected image. The first metric is the relative root-mean-square error (RRMSE), which can be computed by

RRMSE ¼
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi The second metric is the structural similarity (SSIM) index [38], which can be computed by where μ x and μ z are the mean values of imagesx andz, respectively. s 2 x and s 2 z are the variance values of imagesx andz, respectively. σ xz is the correlation coefficient value of imagex and imagez. c 1 = (k1 � s) 2 and c 2 = (k2 � s) 2 where k1 = 0.01 and k2 = 0.03 and s is the dynamic range of the pixel values.

Results of the simulation with noise-free projection data
In this section, we applied the proposed algorithms to the simulated projection data without statistical noise corresponding to the above-mentioned 2D chest slice. The simulated sinogram consisted of 512 (bins) × 1,000 (angles), where the standard parallel-beam geometry with 360å ngular range was assumed. We considered three different cases, each of which corresponds to a different setup of the non-uniform sensitivity bins and consequently a different ring artifact pattern in the image domain. The details are summarized as follows. In the first case, we assumed that some non-adjacent bins have non-uniform sensitivities producing the ring artifact. In the other two cases, some adjacent detector bins have non-uniform sensitivities producing the thick ring artifact of multiple pixel width where the maximum width was set to 3 (pixels). We assumed that the imperfect error value in each detector bin is constant with respect to the angle in the first and second cases, and varied with respect to the angle in the third case. Using the sinogramsb corresponding to the three different cases, we reconstructed imagesx by the proposed algorithms as well as M1 and M2 methods. The two parameters (ρ and δ) control the degree of the ring artifact removal in the proposed algorithms. All the algorithms were implemented using the same values for these parameters ρ, δ, N, and K for the following reasons: to obtain a fair comparison between the different algorithms and to show the effect of the different sparsity-induced norms on the ring artifact elimination. The parameter values used in our implementation are summarized as follows: L = 2, K = 2, 000, α = 125 × 10 −3 , ρ = 0.004, δ = 0.092, b 0 ¼ 15:0, and η = 0.001. Also, the parameter values (σ=0.8, L=4, and W=db40) and (R=5, S m =5, and S l =10) were used for the correction by M1 and M2 methods, respectively. Fig 2 illustrates an example of the reconstructed images using ℓ 2 and ℓ 2smth algorithms for the first study case. The true slice, the corrupted image, the reconstructions using all the proposed algorithms, and the images corrected by using the M1 and M2 methods are shown in Figs 3-5 for all the study cases. Moreover, Figs 6 and 7 show the reconstructed images using the angular constrained forms of the proposed algorithms and their improved versions for the first and second study cases, respectively. Also, an enlargement of a small region of interest (ROI) containing the heart region is shown in the upper right corner of each image with more compressed grayscale. From Fig 2, we can notice that the ring artifact still exists in the reconstructed images using ℓ 2 and ℓ 2 -smth algorithms. On the other hand, it is clear from the reconstructions shown in Figs 3-5 that the M1 method eliminated the ring artifact, but there still exist wave-like errors. Meanwhile, the M2 method partially suppressed the ring artifact. Also, additional artifacts were created like the ones shown in the part highlighted by the red arrows. Choosing different parameters in the M1 and M2 methods may further remove the rings but also damages the image so that we did not succeed in getting nicer images than those shown here. By looking at the reconstructions by the proposed algorithms, we can observe that ℓ 1 -ring algorithm succeeds in eliminating most of the ring artifact except few rings located near the center. For example, we can notice that one ring appeared with low contrast in the first study case, while very small number of rings still appeared in the second and third study cases. Meanwhile, the reconstructions by ℓ 1 -smth-ring algorithm illustrate that the remaining artifacts were smoothed well. The reconstructions by ℓ 0 -ring algorithm produced ring-free images, though very little new artifacts were introduced in the second study case. Additional improvement with the smoothing penalty was able to efficiently suppress these remaining artifacts. It also produced high-quality image closer to the original image without any loss in the fine details. Furthermore, the reconstruction by Hℓ 1 -ring and Hℓ 0 -ring algorithms show that they provide similar results to those in ℓ 1 -ring and ℓ 0 -ring algorithms, respectively. However, the latter algorithms eliminate more ring artifacts and produce better results than the former ones. The residual artifacts were disappeared well with Hℓ 0 -smth-ring algorithm and were smoothed well with Hℓ 1 -smth-ring algorithm, where they still slightly appeared with low contrast. Moreover, it is notable from Figs 6 and 7 that the angular constrained forms of the proposed algorithms were able to suppress more ring artifacts compared to the original non-angular constrained ones. It could eliminate approximately all the ring artifact except very few rings appearing with very low contrast in some algorithms. Further improvement with the smoothing penalty can weaken the remaining artifacts. For quantitative evaluation, we measured the quantitative metrics (RRMSE and SSIM) between the original reference image and the reconstructed image by each algorithm. The obtained metric values are summarized in Tables 2 and 3 for all the study cases. From these results, we can observe that the RRMSE values for the proposed algorithms are significantly smaller than the compared methods. Also, the SSIM values for the proposed algorithms are higher than the compared methods. Furthermore, we note that the metric values for ℓ 1 -smthring and ℓ 0 -smth-ring algorithms and their angular constrained forms outperformed the other proposed algorithms.
Finally, from the above qualitative and quantitative evaluation, the proposed algorithms and their angular constrained forms, especially with incorporating the smoothing penalty, are effective in suppressing the ring artifact. They also achieve higher image quality compared with the other methods implemented in this work. Fig 8 displays both RRMSE and SSIM values versus the iteration numbers for the third study case, and it is clear that the reconstruction

PLOS ONE
by the proposed algorithms required about 2000 iterations to converge. One drawback of the proposed algorithms is that the iterative algorithms derived from ADMM method are slow in convergence. In the future, we plan to improve the slow convergence speed by introducing faster optimization methods.

Results of the simulation with noisy projection data
We also studied the case where projection data contains statistical noise in addition to the ring artifact. In the noisy data case, Poisson noise corresponding to 10 6 (counts/bin) was added to the sinogram of the first study case in the previous section. Using this sinogram, we reconstructed images by the proposed algorithms. The parameter values in the proposed algorithms were as follows: L = 2, K = 500, α = 5 × 10 −3 , ρ = 0.001, δ = 0.03 and η = 0.01. Also, b 0 were 750.0 in ℓ 1 -smth-ring, 55.0 in ℓ 0 -smth-ring, 75.0 in Hℓ 1 -smth-ring, 45.0 in Hℓ 0 -smth-ring, and From the visual evaluation of the results, we can notice that the M1 and M2 methods produced similar results to those in the previously studied noise-free case. Furthermore, we can observe that the reconstructions by ℓ 1 -ring and ℓ 0 -ring algorithms did not succeed in completely eliminating the ring artifact, where it only slightly decreased the intensity of some rings. The reconstructions by Hℓ 1 -ring and Hℓ 0 -ring algorithms produced similar results to those in ℓ 1 -ring and ℓ 0 -ring algorithms, respectively. Meanwhile, the remaining artifacts were weakened with ℓ 1 -smth-ring algorithm. Also, they were suppressed with ℓ 0 -smth-ring, Hℓ 1 - smth-ring, and Hℓ 0 -smth-ring, where the artifact slightly appeared with low contrast. It can be notable that the reconstructions by Hℓ 1 -ring and Hℓ 0 -ring algorithms contained less noise compared to ℓ 1 -ring and ℓ 0 -rings algorithms. The reason for this is the existence of the quadratic part in Huber ℓ 1 -norm and Huber ℓ 0 -norm cost functions in Hℓ 1 -ring and Hℓ 0 -ring algorithms.
Concerning the results of the angular constrained forms of the proposed algorithms, it is clear that the angular constrained forms of the proposed algorithm could eliminate all the ring artifact with ℓ 0 -ring and Hℓ 0 -ring algorithms. They also removed most of the rings except some rings located near the center in both ℓ 1 -ring and Hℓ 1 -ring algorithms, although some new artifacts were introduced. Additionally, we can observe that the angular constrained forms of the improved algorithms could effectively remove the artifact and produce ring-free images. To perform quantitative evaluation of the proposed algorithms, RRMSE and SSIM values were measured and presented in the last two columns in Table 2. We can observe that the improved versions of the proposed algorithms and their angular constrained forms achieved the best values among all the implemented algorithms. The performance of the proposed algorithms in the case of noisy projection is less than that in the case of noise-free projection data, even with Hℓ 1 -ring or Hℓ 0 -ring algorithms. Consequently, we plan to incorporate a statistical model in addition to the sparsity-induced norms for further improvement of the reconstruction with noisy projection data.

Results of the real data
For additional evaluation of the proposed algorithms, we implemented the proposed algorithms with real data. We used the projection data corresponding to the central slice of the above-mentioned real data. We implemented the proposed method and the compared methods with this data. The results are shown in Fig 11. Also, an enlargement of a small ROI located in the center of each reconstruction is shown in the upper right corner of each image. The results are summarized as follows. Most of the ring artifacts were removed well with the corrections using M1 and M2 methods. However, some remaining rings were still visible. These existing techniques worked well, but were not complete in this experiment. The reconstructed images with ℓ 2 and ℓ 2 -smth algorithms and their angular constrained versions (cstr-ℓ 2 and cstr-ℓ 2 -smth) still included the strong ring artifact. Since the real projection data is corrupted with noise, Hℓ 0 -ring, Hℓ 0 -smth-ring, and the angular constraint versions of Hℓ 1 -ring and Hℓ 0 -ring algorithms could remove most of the ring artifact well. Moreover, cstr-smth-Hℓ 1 -ring and cstr-smth-Hℓ 0 -ring algorithms almost completely removed the ring artifact leading to high-quality reconstructions. Finally, we would like to strengthen the following fact. Using ℓ 2 -norm cannot suppress the ring artifact, whereas the combination of ℓ 1 -or ℓ 0 -norm and ℓ 2 -norm as in Huber ℓ 1 -and Huber ℓ 0 -norms is able to remove the ring artifact. This means that the use of sparsity-induced norm is essential in the ring-artifact elimination problem.

Conclusion
In this work, we developed a class of algorithms to eliminate the ring artifact from CT images. The first group of algorithms minimizes different sparsity-induced norms of the imperfect error components of the sinogram. In the second group of algorithms, we improved the first class of algorithms by incorporating the smoothing penalty into the cost function to remove the remaining ring artifact. Furthermore, we also introduced angular constrained forms corresponding to the proposed algorithms and their improved versions. The proposed algorithms are able to correct the ring artifact well during image reconstruction and need neither pre-processing nor post-processing steps. We evaluated the performance of these algorithms through a number of simulation and real data studies. In the simulation, we considered three different cases where different setups of the ring artifact patterns were simulated. Both the qualitative and quantitative evaluation results demonstrated that the improved versions of the proposed algorithms and their angular constrained forms are efficient in ring artifact suppression. They also provide significantly better results than the compared methods, i.e. M1 and M2 methods. In summary, we can conclude that the sparsity-induced norms for the imperfect error components of the projection data in the data fidelity as well as the smoothing penalty can effectively eliminate the ring artifact. Reconstructed images in the case of noisy projection data (noisy projection data, non-adjacent error detector bins, and angular-independent error case).